
# Setup ------------------------------------------------------------------------
  # Options
  options(stringsAsFactors = F)
  # Packages
  library(pacman)
  p_load(data.table, ggplot2, ggthemes, magrittr, extrafont, Cairo)
  # Directories
  dir_project <- "/Users/edwardarubin/Dropbox/Research/MyProjects/NaturalGas/"
  dir_figures <- paste0(dir_project,
    # "Presentations/Beamer20170926HEREW/Figures/")
    "Text/ForSubmission/Figures/")

# Colors and themes ------------------------------------------------------------
  # Colors
  dark_grey    <- "#7D7D7D"
  light_grey   <- "#C7C7C7"
  v_light_grey <- "#EEEEEE"
  mid_grey     <- "#A2A2A2"
  purple       <- "#6A1B9A"
  red_pink     <- "#E91E63"
  orange       <- "#FFA000"
  # Beamer theme
  theme_beamer <- theme_bw() + theme(
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    legend.position = 'bottom',
    legend.key = element_blank(),
    panel.grid = element_blank(),
    text = element_text(family = "Roboto", color = "black",
      size = 16),
    title = element_text(size = 18),
    legend.text = element_text(size = 18),
    legend.key.size = unit(2.5, "cm"))
  # Paper theme
  theme_paper <- theme_bw() + theme(
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    legend.position = 'bottom',
    legend.key = element_blank(),
    panel.grid = element_blank(),
    text = element_text(family = "Charter", color = "black",
      size = 11),
    axis.text = element_text(color = "black"),
    title = element_text(size = 12),
    legend.text = element_text(size = 11),
    legend.key.size = unit(1.5, "cm"))
  # Function for axis decimals
  scale_fun <- function(x) sprintf("%.2f", x)
  # My own save functions
  save_beamer <- function(gg_tmp, name, width = 16 * 0.7, height = 10 * 0.7,
    themes = NULL) {
      ggsave(filename = paste0(name, ".pdf"),
        path = dir_figures,
        plot = gg_tmp + theme_beamer + themes,
        width = width,
        height = height,
        device = cairo_pdf)
      embed_fonts(paste0(dir_figures, name, ".pdf"))
  }
  save_paper <- function(gg_tmp, name, width = 6.4, height = 4,
    themes = NULL) {
      ggsave(filename = paste0(name, ".pdf"),
        path = dir_figures,
        plot = gg_tmp + theme_paper + themes,
        width = width,
        height = height,
        device = cairo_pdf)
      embed_fonts(paste0(dir_figures, name, ".pdf"))
  }

# Supply function --------------------------------------------------------------
  # The function
  supply <- function(x) {
    ifelse(x < 30, 0.9, 1.4)
  }
  # Supply segments (for ggplot2 segments)
  s_dt <- data.table(
    x = c(0, 30, 30), xend = c(30, 30, 50),
    y = c(0.9, 0.9, 1.4), yend = c(0.9, 1.4, 1.4))

# Demand functions -------------------------------------------------------------
  # Summer demand functions (inelastic)
  summer <- function(x) -0.3 * x + 7 + 0.4
  summer_t <- function(x) summer(x) - 0.4
  summer_t2 <- function(x) summer(x) - 0.6
  # Winter demand functions (more elastic)
  winter <- function(x) -0.0025 * x + 0.95 + 0.4
  winter_t <- function(x) winter(x) - 0.4
  winter_t2 <- function(x) winter(x) - 0.2

# Solve for points: Summer -----------------------------------------------------
  # Find intersection between first tier and summer_t
  s1 <- data.table(
    x = uniroot(function(x) summer_t(x) - 0.9, c(-1e3,1e3))$root,
    y = 0.9)
  # Find intersection between first tier and summer_t2
  s1_2 <- data.table(
    x = uniroot(function(x) summer_t2(x) - 0.9, c(-1e3,1e3))$root,
    y = 0.9)
  # Find intersection between first tier and summer
  s2 <- data.table(
    x = uniroot(function(x) summer(x) - 0.9, c(-1e3,1e3))$root,
    y = 0.9)
  # If summer demand hits the step, then calculate the point of intersection
  if (between(summer(30), 0.9, 1.4)) {
    s3 <- data.table(
      x = 30,
      y = summer(30))
    } else {
      s3 <- NULL
    }
  # If summer_t hits the step, then calculate the point of intersection
  if (between(summer_t(30), 0.9, 1.4)) {
    s4 <- data.table(
      x = 30,
      y = summer_t(30))
    } else {
      s4 <- NULL
    }
  # Calculate the vertical projection of taxed demand on untaxed demand
  s5 <- data.table(x = s1$x, y = summer(s1$x))
  s5_2 <- data.table(x = s1_2$x, y = summer(s1_2$x))
  # Update s2$x, if s3$y is at the step between tiers
  if (!is.null(s3)) if (between(s3$y, 0.9, 1.4)) s2$x <- 30
  # Join the points
  summer_dt <- rbindlist(list(s1, s2, s3, s4, s5))
  summer_dt[, id := as.factor(1)]
  summer_dt2 <- rbindlist(list(s1_2, s5_2, s5, s1))
  summer_dt2[, id := as.factor(2)]

# Solve for points: Winter -----------------------------------------------------
  # Find intersection between first tier and winter_t
  w1 <- data.table(
    x = uniroot(function(x) winter_t(x) - 0.9, c(-1e3,1e3))$root,
    y = 0.9)
  # Find intersection between first tier and winter
  w2 <- data.table(
    x = uniroot(function(x) winter(x) - 0.9, c(-1e3,1e3))$root,
    y = 0.9)
  # If winter demand hits the step, then calculate the point of intersection
  if (between(winter(30), 0.9, 1.4)) {
    w3 <- data.table(
      x = 30,
      y = winter(30))
    } else {
      w3 <- NULL
    }
  # If winter_t hits the step, then calculate the point of intersection
  if (between(winter_t(30), 0.9, 1.4)) {
    w4 <- data.table(
      x = 30,
      y = winter_t(30))
    } else {
      w4 <- NULL
    }
  # Calculate the vertical projection of taxed demand on untaxed demand
  w5 <- data.table(x = w1$x, y = winter(w1$x))
  # Update w2$x, if w3$y is at the step between tiers
  if (!is.null(w3)) if (between(w3$y, 0.9, 1.4)) w2$x <- 30
  # Join the points
  winter_dt <- rbindlist(list(w1, w2, w3, w4, w5))
  winter_dt[, id := as.factor(1)]
# NOTE: These calculations assuming the new tax line hits the step
  # Project winter_t quantity onto winter_t2
  w2_b <- data.table(x = w1$x, y = winter_t2(w1$x))
  # Find where winter_t2 hits the step
  w3_b <- data.table(x = 30, y = winter_t2(30))
  # The corner
  w4_b <- data.table(x = 30, y = 0.9)
  winter_dt2 <- rbindlist(list(w1, w2_b, w3_b, w4_b))
  winter_dt2[, id := as.factor(2)]

# Plot summer: current regime --------------------------------------------------
  gg_tmp <- ggplot() +
    # DWL triangle
    geom_polygon(data = summer_dt, aes(x, y, group = id), fill = light_grey) +
    # Demand: CARE in summer, without tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = summer) %>% unlist()),
      aes(x = x, y = y), color = mid_grey, linetype = 2) +
    # Demand: CARE in summer, with tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = summer_t) %>% unlist()),
      aes(x = x, y = y)) +
    # Supply
    geom_segment(data = s_dt, aes(x = x, xend = xend, y = y, yend = yend)) +
    geom_point(data = s1, aes(x,y), size = 2) +
    xlim(c(0,50)) + xlab("") +
    scale_y_continuous(
      name = "",
      limits = c(0.7,1.6),
      labels = scale_fun) +
    # Axes
    # geom_vline(xintercept = 0, color = dark_grey, size = 0.4) +
    geom_vline(xintercept = 0, color = "black", size = 0.6) +
    # Label lines
    annotate(geom = "text", x = 49.5, y = 1.43,
      # label = "S", size = 7.5, family = "Robot") +
      label = "S", size = 4, family = "Georgia") +
    annotate(geom = "text", x = 16.6, y = 0.75,
      # label = "D[tau]", size = 7.5, parse = T, family = "Roboto")
      label = "D[tau]", size = 4, parse = T, family = "Georgia")
  # Save
  # save_beamer(gg_tmp, "dwlSummerNow", width = 5, height = 6.25)
  save_paper(gg_tmp, "dwlSummerNow", width = 3, height = 4)

# Plot summer: new regime ------------------------------------------------------
  gg_tmp <- ggplot() +
    # DWL triangle
    geom_polygon(data = summer_dt, aes(x, y, group = id), fill = light_grey) +
    geom_polygon(data = summer_dt2, aes(x, y, group = id), fill = red_pink) +
    # Demand: CARE in summer, without tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = summer) %>% unlist()),
      aes(x = x, y = y), color = mid_grey, linetype = 2) +
    # Demand: CARE in summer, with tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = summer_t) %>% unlist()),
      aes(x = x, y = y), color = light_grey, alpha = 0.5) +
    # Demand: CARE in summer, with new tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = summer_t2) %>% unlist()),
      aes(x = x, y = y)) +
    # Supply
    geom_segment(data = s_dt, aes(x = x, xend = xend, y = y, yend = yend)) +
    geom_point(data = s1_2, aes(x,y), size = 2) +
    xlim(c(0,50)) + xlab("") +
    scale_y_continuous(
      name = "",
      limits = c(0.7,1.6),
      labels = scale_fun) +
    # Axes
    # geom_vline(xintercept = 0, color = dark_grey, size = 0.4) +
    geom_vline(xintercept = 0, color = "black", size = 0.6) +
    # Label lines
    annotate(geom = "text", x = 49.5, y = 1.43,
      # label = "S", size = 7.5, family = "Roboto") +
      label = "S", size = 4, family = "Georgia") +
    annotate(geom = "text", x = 16.1, y = 0.75,
      # label = "D[tau*minute]", size = 7.5, parse = T, family = "Roboto")
      label = "D[tau*minute]", size = 4, parse = T, family = "Georgia")
  # Save
  # save_beamer(gg_tmp, "dwlSummerNew", width = 5, height = 6.25)
  save_paper(gg_tmp, "dwlSummerNew", width = 3, height = 4)

# Plot winter: current regime --------------------------------------------------
  gg_tmp <- ggplot() +
    # DWL triangle
    geom_polygon(data = winter_dt, aes(x, y, group = id), fill = light_grey) +
    # Demand: CARE in winter, without tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = winter) %>% unlist()),
      aes(x = x, y = y), color = mid_grey, linetype = 2) +
    # Demand: CARE in winter, with tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = winter_t) %>% unlist()),
      aes(x = x, y = y)) +
    # Supply
    geom_segment(data = s_dt, aes(x = x, xend = xend, y = y, yend = yend)) +
    geom_point(data = w1, aes(x,y), size = 2) +
    xlim(c(0,50)) + xlab("") +
    scale_y_continuous(
      name = "",
      limits = c(0.7,1.6),
      labels = scale_fun) +
    # Axes
    # geom_vline(xintercept = 0, color = dark_grey, size = 0.4) +
    geom_vline(xintercept = 0, color = "black", size = 0.6) +
    # Label lines
    annotate(geom = "text", x = 49.5, y = 1.43,
      # label = "S", size = 7.5, family = "Roboto") +
      label = "S", size = 4, family = "Georgia") +
    annotate(geom = "text", x = 49.1, y = 0.86,
      # label = "D[tau]", size = 7.5, parse = T, family = "Roboto")
      label = "D[tau]", size = 4, parse = T, family = "Georgia")
  # Save
  # save_beamer(gg_tmp, "dwlWinterNow", width = 5, height = 6.25)
  save_paper(gg_tmp, "dwlWinterNow", width = 3, height = 4)

# Plot winter: new regime ------------------------------------------------------
  gg_tmp <- ggplot() +
    # DWL triangle
    geom_polygon(data = winter_dt, aes(x, y, group = id), fill = light_grey) +
    geom_polygon(data = winter_dt2, aes(x, y, group = id), fill = purple) +
    # Demand: CARE in winter, without tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = winter) %>% unlist()),
      aes(x = x, y = y), color = mid_grey, linetype = 2) +
    # Demand: CARE in winter, with tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = winter_t) %>% unlist()),
      aes(x = x, y = y), color = light_grey, alpha = 0.5) +
    # Demand: CARE in winter, with new tax
    geom_line(data = data.table(
      x = seq(0, 50, 0.1),
      y = lapply(X = seq(0, 50, 0.1), FUN = winter_t2) %>% unlist()),
      aes(x = x, y = y)) +
    # Supply
    geom_segment(data = s_dt, aes(x = x, xend = xend, y = y, yend = yend)) +
    geom_point(data = w2_b, aes(x,y), size = 2) +
    xlim(c(0,50)) + xlab("") +
    scale_y_continuous(
      name = "",
      limits = c(0.7,1.6),
      labels = scale_fun) +
    # Axes
    # geom_vline(xintercept = 0, color = dark_grey, size = 0.4) +
    geom_vline(xintercept = 0, color = "black", size = 0.6) +
    # Label lines
    annotate(geom = "text", x = 49.5, y = 1.43,
      # label = "S", size = 7.5, family = "Roboto") +
      label = "S", size = 4, family = "Georgia") +
    annotate(geom = "text", x = 49.1, y = 1.06,
      # label = "D[tau*minute]", size = 7.5, parse = T, family = "Roboto")
      label = "D[tau*minute]", size = 4, parse = T, family = "Georgia")
  # Save
  # save_beamer(gg_tmp, "dwlWinterNew", width = 5, height = 6.25)
  save_paper(gg_tmp, "dwlWinterNew", width = 3, height = 4)
